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^ : ABSTRACT 
O 

^ I We calculate numerically the collapse of slowly rotating, non- magnetic, mas- 

O I sive molecular clumps of masses 30 Mq, 60 Mq, and 120 M©, which conceivably 

could lead to the formation of massive stars. Because radiative acceleration on 
dust grains plays a critical role in the clump's dynamical evolution, we have 



S ■ improved the module for continuum radiation transfer in an existing 2D (axial 



6 ■ 

^ ■ symmetry assumed) radiation hydrodynamic code. In particular, rather than 

using "grey" dust opacities and "grey" radiation transfer, we calculate the dust's 
^ wavelength-dependent absorption and emission simultaneously with the radia- 

' tion density at each wavelength and the equilibrium temperatures of three grain 

components: amorphous carbon particles, silicates, and "dirty ice" -coated sili- 
cates. Because our simulations cannot spatially resolve the innermost regions 
of the molecular clump, however, we cannot distinguish between the formation 
of a dense central cluster or a single massive object. Furthermore, we cannot 
exclude significant mass loss from the central object (s) which may interact with 
the inflow into the central grid cell. Thus, with our basic assumption that all 
material in the innermost grid cell accretes onto a single object, we are only able 
to provide an upper limit to the mass of stars which could possibly be formed. 

We introduce a semi-analytical scheme for augmenting existing evolution- 
ary tracks of pre-main sequence protostars by including the effects of accretion. 
By considering an open outermost boundary, an arbitrary amount of material 
could, in principal, be accreted onto this central star. However, for the three 
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cases considered (30 Mq, 60 Mq, and 120 Mq originally within the computa- 
tion grid), radiation acceleration limited the final masses to 31.6 Mq, 33.6 Mq, 
and 42.9 Mq, respectively, for wavelength-dependent radiation transfer and to 
19.1 Mq, 20.1 Mq, and 22.9 Mq for the corresponding simulations with grey 
radiation transfer. Our calculations demonstrate that massive stars can in prin- 
ciple be formed via accretion through a disk. The accretion rate onto the central 
source increases rapidly after one initial free-fall time and decreases monotoni- 
cally afterwards. By enhancing the non-isotropic character of the radiation field 
the accretion disk reduces the effects of radiative acceleration in the radial direc- 
tion — a process we denote the "fiashlight effect" . The fiashhght effect is further 
amplified in our case by including the effects of frequency dependent radiation 
transfer. We conclude with the warning that a careful treatment of radiation 
transfer is a mandatory requirement for realistic simulations of the formation of 
massive stars. 

Subject headings: hydrodynamics — radiation transfer — stars: formation — 
circumstellar disks — outflows 



1. Introduction 

Although massive stars play a critical role in the production of turbulent energy in the 
ISM, in the formation and destruction of molecular clouds, and ultimately in the dynami- 
cal and chemo-dynamical evolution of galaxies, our understanding of the sequence of events 
which leads to their formation is still rather limited. Because of their high luminosities we 
can expect: a) radiative acceleration will contribute significantly to the dynamical evolution 
during the formation process and b) the thermal evolution time scales of massive pre-main 
sequence objects will be extremely short. Thus, we cannot simply "scale up" theories of low 
mass star formation. Furthermore, OB stars form in clusters and associations; their mu- 
tual interactions via gravitational torques, powerful winds and ionizing radiation contribute 
further to the complexity of the problem. 

Even though no massive disk has yet been directly observed around a main sequence 
massive star, it is likely that such disks arc the natural consequence of the star formation 
process even in the high mass case. In their radio recombination maser studies and CO 
measurements Martin-Pintado et al. (1994) do find indirect evidence for both an ionized 
stellar wind and a neutral disk around MWC349. Moreover, several other high luminosity 
FIR sources — suspected embedded young OB stars — have powerful bipolar outfiows 
associated with them (e.g., Eiroa et al. 1994; Shepherd et al. 2000). Such massive outfiows 
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are probably powered by disk accretion, and, similar to their low mass counterparts, the 
flow energetics appear to scale with the luminosity of the source (see Cabrit & Bertout 1992; 
Shepherd & Churchwell 1996; Richer et al. 2000). 

The detailed structure and evolutionary history of massive circumstellar disks has im- 
portant consequences with regard to the early evolution of these protostars. Disks provide 

a reservoir of material with specific angular momentum too large to be directly accreted by 
the central object. Only after angular momentum is transported outwards can this material 
contribute to the final mass of the star. The transition region disk-star will strongly infiuence 
the star's photosphcric appearance and how the star interacts with the disk. The relative 
high densities in these disks provide the environment for the further growth and evolution of 
dust grains, affecting the disk's opacity and consequently its energetics and appearance. The 
disk can be expected to interact with stellar outfiows and is likely to be directly responsible 
for the outfiows associated with star formation. 

Disks surrounding massive stars or disks associated with close companions to massive 
stars should be short-lived compared to their low mass counterparts. The UV environment 
within an OB cluster will lead to the photoevaporation of disks on a time scale of a few 
10^ yr (HoUenbach, Yorke, & Johnstone 2000). Because this process operates on a time scale 
comparable to the formation of massive stars and is competitive to it, it is important to 
carefully model the transfer of radiation in the envelopes of accreting massive stars. Numer- 
ical tools capable of this task are lacking at present. We consider the present investigation 
as an important step in this direction. 

2. The numerical model 

We consider the simulation of the hydrodynamic collapse of a rotating molecular cloud 
clump, with wavelength-dependent radiation transfer, under the assumption of symmetry 
with respect to the rotation axis and the equatorial plane. Our code contains all of the basic 
features of the 2-D code described in detail by Yorke & Bodenheimer (1999; hereafter YB). 
In the following we shall discuss only the deviations from and enhancements to the YB code. 

In addition to artificial viscosity for the treatment of shocks, physical viscosity has been 
implemented via an a prescription (Shakura & Sunyaev 1973): v = acs{r)H{r), where 
we have approximated the local disk scale height by H{r) = cs{r)/fl{r) {fl is the angular 
velocity and cs is the sound speed). All components of the viscosity tensor are included, 
as previously implemented in 2-D disk models by Rozyczka et al. (1994) and YB. Contrary 
to YB, who allowed a to vary in time in order to mimic the effects of tidal torques due to 
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the growth of non-axisymmetric gravitational modes, the value for a = 10 ^ is was kept 
constant in space and time for all cases considered. 

The details of this angular momentum transfer scheme — within certain limits — do 
not critically affect our results, however, because our central computational zone (where 
angular momentum transfer is presumably very critical) is so large. Test calculations with 
a = 0.03 yielded essentially the same results. By contrast, for a < 10~^ ring instabilities 
developed in the disk which led to a premature ending of the calculations. As discussed by 
Yorke, Bodenheimer, Laughlin (1995), who did not include the effects of angular momentum 
transfer (i.e., a = 0), these ring instabilities arc unrealistic. Such a ring would be unstable 
on a very short (dynamical) time scale and the resulting clumps would exert tidal torques 
resulting in angular momentum transfer. 

As in YB we utilize a series of hierarchically nested grids (Berger & Colclla 1989; Yorke 
& Kaisig 1995). Whereas YB considered nesting levels of 6 and 60 x 60 grids, allowing the 
innermost grid to be ~ 1/2000 of the cloud radius, here we have considered nesting levels 
of 3 only and slightly larger grids (64 x 64). We assume constant density p = po along the 
outermost radius Tmax and allow material to enter into or exit from our computational grid 
R'^ + Z'^ < r^g^ (R and Z are the cyhndrical coordinates), based on the sign of the radial 
component of the velocity. By contrast, YB assume a semi-permeable outer boundary at 
Tinax^ Material with positive radial velocity can leave the computational grid but no mass 
was allowed to enter. 

As in YB the boundary of the the innermost cell of our innermost grid is considered 
to be semi-permeable: Material can flow into this cell but cannot flow out of it. Material 
entering this cell is assumed to accrete onto a single central object. We realize that this 
"sink cell" procedure is a gross simplification of the physics in the innermost regions of 
our computational domain and a number of possibly important effects are being ignored, 
e.g. fragmentation and accretion of material onto multiple objects and the interaction of 
the accretion flow with powerful outflows. For the cases F30 and G30 (see section 3), for 
instance, our innermost cell is a cyhnder of radius 40 AU and height 80 AU, whereas for the 
cases F120 and G120 the cell is larger by a factor of four. Although we cannot follow the 
mass flow within and possibly out of this cell, our simulations do provide upper limits to 
the amount of material which is available to be accreted by the central object: If radiative 
acceleration prevents the flow of material into the central sink cell, then the central object 
cannot accrete it. 



- 5 - 



2.1. Modeling the Central Star 

In contrast to YB we consider a slightly more sophisticated treatment of the radius i?*, 
luminosity L*, and effective temperature Tcs of the central object. Its mass is uniquely 
determined by integrating the mass flux into the center sink cell: 

= j M^dt. (1) 

We can express the total energy of the central core in terms of a "structure parameter" tj: 

^« = -V^ . (2) 

In principal, 77, a parameter describing the compactness of the hydrostatic core, must be 
calculated by solving for the stellar structure of the accreting hydrostatic core. For polytropes 
of degree n, rj can be derived analytically (e.g. Kippenhahn & Weigert 1990): 

A fully convective pre-main sequence protostar can be approximated by an n = 3/2 polytrope 
and 1] = 3/7. As the star approaches the main sequence, a greater proportion of it becomes 
radiative, its core becomes more compact, and r] increases. 

For purposes of discussion we shall assume for the moment that 77 = r]{M^^R^) is a 
known function. In this case the intrinsic core luminosity (which does not include the 
contribution to the total luminosity emitted in the accretion shock front and dissipated 
within the disk) is given by: 

L* — i^nuc ~ -C'tot ~ y' 



R 

-^nuc ~ -E'tot 



!Z + f 2 - ^ 1 ^ - ^ 

7] \ r]jM^R^ 



(4) 



where Lnuc is the contribution from nuclear burning. 



In the following we will assume that P — 1: The material accreted by the star is being 
added ever so gently at its current radius i?* and that this material adds negligible entropy 
to the star.^ For the total bolometric luminosity of a spherically accreting star, however, we 



"'^For the newly accreted material we nmst subtract the difference of potential energy from infinity to the 
stellar surface when considering the star's net change of total energy. For /3 = 1 the energy gained by the 
star due to heating from the accretion shock ( "backwarming" ) or by dissipating rotational energy within the 
star is negligible. 
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must include the contribution of the potential energy of infalling material as it is dissipated 
on its way to the stellar surface: 

T - T T - T nGM.M, 

-^bol — -L/* + -i^acc — -Li* + P — (Oj 

-ft* 

Equation 4 shows that for constant — Lnuc the star should increase or decrease its 
radius due to mass accretion, depending on the sign of the coefficient of M^,/M^,. Because 
(3^1 and rj ^ 3/7 during the fully convective Hayashi phase, this coefficient is negative 
(~ —1/3). Thus, mass accretion onto a fully convective star has a tendency to decrease the 
star's radius, whereas close to the main sequence, where rj is closer to unity, mass accretion 
has a tendency to cause the star to bloat up (see Kippenhahn & Hofmeister 1977 for a more 
detailed discussion). In reality, the internal readjustment of the star after it has gained mass 
also affects the nuclear burning rate and thus has an effect on both radius and luminosity, 
depending on the magnitude of the mass accretion rate and the star's current position in the 
Hertzsprung-Russell (HR) diagram. 

A classical problem of the mathematical theory of stellar structure is the question of 
whether for stars of given fixed parameters, say chemical composition, mass, and radius, 
there exists one and only one solution of the basic structure equations. There is actually 
no mathematical basis for the so-called "Vogt- Russell" conjecture of uniqueness and indeed, 
multiple solutions for the same set of parameters have been found numerically in some 
cases (see discussion by Kippenhahn & Wcigcrt 1990). However, for the rather simple cases 
considered here, spherically symmetric, quasi-hydrostatic homogeneous prc-main sequence 
and young main sequence stars, knowledge of the star's mass and age at any given time does 
allow us to uniquely fix its position in the HR diagram, from which we can determine L*, 

Teff, and -Lnuc- 

For the pre-main sequence phase we shall account for deuterium burning only and use 
the following approximate expression: 

Lnnc ^ Lo{M^) 

where Lq and Rq are the equilibrium deuterium burning rate and equilibrium radius for a 
star of mass at its "birthline",^ Xofl is the cosmic mass abundance of deuterium, and 
Xd is the star's net deuterium abundance. We have assumed p = 21, which — because the 



^ There are alternate definitions of the concept of "birthhne". Here we use the word to describe the 
equiUbrium position of a homogeneous, deuterium-burning, pre-main sequence star with a cosmic abundance 
of deuterium. 



Xd 

XDfi 



R 



(6) 
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star's central density pc oc R ^ — corresponds to Ld oc p^. This insures that a non-accreting 
star remains close to its birthline until a significant fraction of its deuterium is consumed. 

Assuming instantaneous mixing during accretion, the deuterium mass fraction xd can 
be calculated from the following equation: 

— — = Xd,oM^ - €dLd , (7) 

where eoLi) is the rate of deuterium consumption due to deuterium burning [sd = 1-76 x 
10"^^s^cm~^ is a constant). 

Prom equations 4 and 6 we can derive an expression for R^: 
R* 



JrC^ 

2 - - + 7]M 



M^: — Lnuc 1 /o\ 



tot 



where 



Prom knowledge of ri{M^,R^), L^{M^,R^), Lq{M^), and Rq{M^) we approximate the pre- 
main sequence evolution of an accreting protostar by integrating equations 7 and 8 simulta- 
neously. 

How does one actually determine 77 and from knowledge of M* and i?*? Por the 
Hayashi phase we have used published pre- main sequence tracks^ for L*(M*,r), R^{M^,t) 
and set rj = 3/7. Por the mass range 0.1 Mq < M* < 2.5 M0 we use the evolutionary tracks 
of D'Antona & Mazzitelli (1994), which assume "CM convection" (q;ml = 2) and Alexander 
+ RI opacities (F = 0.28, Z = 0.019). For 3 < M, < 15 Mq we use tracks pubhshed by 
Iben (1965). Both sets of tracks represent a series of stellar models which incorporate the 
detailed microphysics of convection and stellar atmospheres. Por masses M^, > 15 M© we 
have assumed that the evolutionary tracks for non-accreting stars are horizontal lines at the 
main sequence luminosity, similar to the 15 Mq track (see Pig. 1), and the main sequence 
values of R^{M^) and L*(M*) were taken from Allen (1973). 

As the contracting protostellar becomes radiative, 77 increases. Once the star reaches the 
main sequence hydrogen burning commences and 77 ~ const. We approximate 77 during this 



^We use the evolutionary time r given for the published tracks as a parameterization of the curves. Our 
evolutionary time t results from integrating equation 8. 
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phase by approximating the main sequence by polytropic models with the restriction that the 
main sequence value rjus ^ 1- Obviously, this is a very rough approximation, but because 
we expect rj ^ and Lnuc on the main sequence, we are not making a significant error. 
Our greatest error arises during the transition from 77 = 3/7 (fully convective) to rjMS, where 
we have used interpolated values. 

Of course, we could obtain a much better approximation to these physical parameters by 
solving the full set of radiation hydrodynamic equations for an accreting hydrostatic object. 
This, however, goes far beyond the scope of the present investigation. We merely wish to 
obtain more realistic approximations for and i?* than those used by YB which were based 
on the pre-niain sequence evolution of non-accreting low mass protostars. Because massive 
stars require at least one phase of high accretion rates, typically > 10~^ Mq yr"^, 
the effects of accretion on the evolution of and cannot be completely neglected. To 
illustrate this quantitatively we display the evolution of (proto-) stars accreting at given 
constant rates in the HR diagram (Fig. 1). 

These tracks compare very well with published, more detailed calculations by Behrend 
& Maeder (2001) and by Meynet & Maeder (2000). Not only do the tracks lie shghtly below 
the equilibrium deuterium burning "birthline" in all cases, but the qualitative effect of rapid 

accretion — namely to shift the tracks to even smaller radii below the "birthline" — occurs 
both in our simplified model and the above cited calculations. The reason for this is the 
negative sign of 2— (3/7] for fully convective stars. The tracks of our simplified model converge 
to the main sequence in a manner similar to those of the detailed calculations. 

There are some differences, however. Behrend & Maeder and Meynet & Maeder begin 
their tracks assuming a fully convective 0.7 M© star taken to be 7 x 10^ yr old, whereas we 
show here that the equilibrium deuterium-burning position at cosmic deuterium abundance 
is never reached via accretion for masses M < 1 Mq. If you add mass too quickly, the star's 
radius is reduced as discussed above. If you add it too slowly (see 10^^ M© yr~^ track in 
Fig. 1), a significant amount of the deuterium is consumed even before 0.7 M© has been 
accreted. Other differences are accountable by the different time dependent accretion rates. 

We remind the reader that these tracks in the HR diagram do not refiect the actual 
observable bolometric luminosities of accreting protostars. Much of the accretion luminosity 
will be indistinguishable from the intrinsic luminosity of the star. For our hydrodynamic 
simulations we will include the effects of the accretion luminosity when discussing the star's 
evolution within the HR diagram. 

As in YB we add the accretion luminosity Lace (Adams & Shu 1986) to the core's 
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Fig. 1. — Evolutionary tracks of pre- main sequence stars accreting material at a constant 
rate as indicated. All accreting tracks begin at M^, = 0.1 Mq and i?* = 1 R©. The total 
bolometric luminosity, which includes contributions Lace from the relaxation zone behind 
the accretion shock and emission from an accretion disk, is not shown. For comparison, the 
evolutionary tracks of non-accreting stars according to D'Antona & Mazzitelli (1994) and 
Iben (1965) are given {solid grey lines); the tracks for 0.2 Mq, 1.0 Mq, 3.0 Mq, and 15 Mq 
are drawn slightly thicker. 
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intrinsic luminosity to obtain the total luminosity 




3 GM,M, 

4 R, 



(10) 



Equation 10 differs from equation 5 (for f3 = 1), because approximately 1/4 of the total po- 
tential energy of the accreted material is dissipated within the disk and is already accounted 
for by our treatment of viscosity. 

Finally, knowing Ltot and allows us to determine T^e from: 



(TsB is the Stefan-Boltzmann radiation constant. When discussing the evolution of the central 
star in the HR diagram, we use L* rather than Ltot in equation 11. 



The principal source of opacity is due to absorption and scattering by dust (cf. Yorke & 
Henning 1994). We have adopted the detailed frequency dependent grain model of Preibisch 
et al. (1993), which assumes a mixture of small amorphous carbon particles (for grain 
temperatures T^c < 2000 K) and "astrophysical silicate" grains (for temperatures Tsi < 
1500 K; cf. Draine & Lee 1984). At grain temperatures Tsa < 125 K the silicates are coated 
with a layer of "dirty" NH3/H2O ice, contaminated with 10% of the amorphous carbon 
particles. The grain sizes are assumed to follow an MRN power law n(a) ~ a~^'^ (Mathis, 
Rumpl, & Nordsieck, 1977) in the size ranges 7nm < Cac < SOnm (carbon particles) and 
40nm < osi < l/^m (silicates). As evident in Figure 2, the specific extinction is strongly 
wavelength dependent. 

For all our calculations we used 64 non-uniformly distributed frequency points between 
5000 /im and 0.1 /im. 



A complete description of the radiation field requires knowledge of the radiation inten- 
sity /jy(x, n, t), where x indicates the position and fi the direction under consideration. Even 
assuming an axially symmetric configuration and plane symmetry with respect to the equa- 
torial plane, I^, is a function of 6 independent variables: frequency i/, two spatial variables 



'tot — 



ATiasBR'Xs ■ 



(11) 



2.2. The Opacity Model 



2.3. Modeling Radiation Transfer 
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Wavelength [ jim ] 




Frequency [ Hz 



Fig. 2. — Specific extinction a'^^ of tlie dust at temperatures below 125 K (aC + Silicate + 
Ice), between 125 K and 1500 K {aC + Silicate), and between 1500 K and 2000 K (aC). 
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{R, Z), two direction variables, say {9, 0), and time t. Solving the rather innocuous looking 
equation for radiation transfer 

^^ + V-/. = CM'5.-/.) (12) 

{Su is the source function and k^^* = J2i i^Ti is the net contribution to the extinction co- 
efficient from all components "i") together with the equations for hydrodynamics, energy 
balance, radiation equilibrium, and the Poisson equation for the gravitational potential be- 
comes an formidable task with present day computers due to the vast number of computations 
which are necessary to obtain reasonable numerical resolution. 



2.3.1. Flux-limited diffusion 

The numerical problem can be greatly simplified by resorting to the fiux-limited diffusion 
(FLD) approximation (see Yorke & Kaisig 1995, Lcvcrmore & Pomraning 1981). If we denote 
by Jv and H,^ the zeroth and first moments of the radiation field, respectively, where 

Jv—-^ ludO. and Hj, — ^ / I^hdQ, , (13) 

Att J An J 

then the time independent form of the zeroth moment equation — obtained by integrating 
equation 12 over all directions — can be written: 

V • = J, + 5^ [B,{T,) - J,] . (14) 

i 

K^l is the absorption coefficient of grain type "i", Tj is its temperature, By{T) is the Planck 
function, is the contribution to the emissivity from internal sources other than dust emis- 
sion. The next higher order moment equation, obtained by multiplying the time independent 
form of equation 12 by n and then integrating over all directions, would contain terms in- 
volving the second moment of the radiation field. 

The FLD approximation is a procedure for closing this series of moment equations 
by approximating the second moment in terms of the zeroth and first moments, utilizing 
knowledge of the opacity distribution: 

^ -^^-^^ ■ ^'^^ 

The effective albedo (v^, and the flux limiter A^, are deflned as (/t^'Jf is the scattering coefficient 
for grain type "i" ) : 

i^abs D {rp X _i_ sea j 



-13- 



= i ( cothd - -/\ with d = ^^/"^ (17) 

In the hmits of high opacity or low opacity, the FLD approximation asymptotically ap- 
proaches the diffusion limit or the free-streaming limit, respectively, as expected. 

The equations 14 and 15 have to be solved for all frequencies simultaneously with the 
condition for radiative equilibrium of each of the absorbing species: 



(18) 



When considering "grey" radiation transfer as YB did, this condition can be put into tabular 
form (for rapid look-up): 

Ti = %{J) with J = J J^du . (19) 



2.3.2. Finite Difference Equations 

Following a "staggered mesh" discretization philosophy for converting partial derivatives 
into finite differences, we define J^, at the cell centers of an {Rj, Z^) cylindrical grid and the 
vector — {Hj^^r, Hj,^z) at the centers of cell boundaries. Combining the equations 14 and 
15 we obtain a diffusion equation 

V • (p.v J.) = J. + ^ </ [J. - s,(r,)] . (20) 

i 

The discretized form of this FLD equation can be expressed as a matrix equation: 

AJ. = e.(T) , (21) 

where the vector J,^ represents the solution for the zeroth moment and e,^(T) the dust 
temperature dependent components for the source terms at all grid centers. We employ an 
equidistant (i?j, Z^) grid at each level of nesting and, while utilizing a 5-point discretization 
scheme, 

+(^'lJ,kJl,j+l,k + CI'l'j,kJl,j,k+l 

-ji,j,k-Ei>^t^j,kBimj,k), (22) 

insure by proper centering that our finite difference equation 21 is accurate to second order 
0{h^) of the grid spacing h — Rj+i — Rj — Zk+i — Zk- In equation 22 we have used the 
subscript I for frequency and i for grain type. 
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Note that Ai^ is an implicit function of J^. Each element of Au for frequency at a 
grid cell (j, k) contains, with proper centering, J,^/ dependencies for all frequencies z/' and 
extending beyond the 5-point discretization star of grid cells {j,k) and {j ± ib 1). As 
evidenced by equation 18 the right hand side of equation 22 is also an implicit function 
of Jjy. The fact that we must find a new solution iteratively on each nested grid for each 
hydrodynamic time step places strong demands on our solution algorithm. 



2.3.3. Boundary conditions 

Boundary conditions are required for each level of nested grids. For each grid cell 
{Rj,Zk) on the outermost grid which satisfies the condition R'j + > r^^^^^ we specify 
Jjy = w B,j{Tout), where Tout is the color temperature of an external isotropic radiation field 
and w is the radiation dilution factor. For all cases considered here we choose w = 1 and 
Tout = 20 K. Along the outer edges of interior (fine) grids we use the interpolated value of J„ 
from the next (coarser) grid level. Additional boundary conditions along the rotation axis 
and equatorial plane result from the assumed symmetry. 

Within the innermost central grid cell of each level of nested grids we treat the central 
star as an additional internal source of emissivity and define ji, (see equation 14) accordingly: 



= TTZ^MTetl) , (23) 



27rRlZi 

where Ri and Zi are the radial extent and height of the central cell. 



2.4. Solution algorithms 

Because of the necessity of solving equations 18 and 20 repeatedly on several grids 
during the course of hydrodynamic evolution, it was imperative to make this module fast 
and robust. Several promising algorithms which rely on fine-tuning adjustable parameters 
had to be abandoned, because a wide range of problem classes (optically thin, optically thick, 
strong density gradients, emission-dominated, scattering-dominated, etc.) occurred, which 
were not well suited to a single set of parameters. We were able to make vast improvements 
with respect to the iterative procedures used by Sonnhalter et al. (1995), who solved these 
same frequency dependent equations for a few selected density configurations. We sometimes 
sacrificed robustness (but never accuracy) for speed but automatically fell back to slower and 
more robust iterative schemes when the "faster" algorithms failed to converge. We feel it 
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is useful to discuss our "failures" as well as our "successes" in our endeavor to improve the 
speed of our frequency dependent radiation transfer module. 

2.4- 1- Temperature determination 

In analogy to the transfer of line radiation in stellar atmospheres we consider the method 
of approximate A operators (see e.g. Cannon 1973a,b and Scharmer 1981). Rewriting equa- 
tion 21 with the operator A = A~^, we find the formal solution: 

J, = A,e,(T) (24) 

A simple A-iteration would entail calculating T"^'^ and thus ei,(T°''^) from J"''^ using equation 
18. Prom this, a new, improved estimate for J^, can be calculated: J"^^ = A,^e,^(T°''^). 
Replacing J°^'^ by J"*^^ and repeating this procedure several times may or may not (usually 
not) quickly converge to the equilibrium values for T and Ji, satisfying equation 24. This 
procedure is labelled "A" in Table 1. 

In order to improve convergence we consider an appropriate approximation A* to our 
operator A,^. Rather than using J"'*^ in equation 18 we substitute the expression 

jnew ^ joid ^ [e^(T"^") - e,(T°i<^)] (25) 
to obtain the equilibrium conditions for each grain component "i" : 

= /</[j°i^-A:e.(T°i<i)] dv, (26) 

which can be rewritten in the form G(T°'^^) = 0. The solution vector T"^^ can be obtained 
by solving equation 26 using a multidimensional Newton-Raphson procedure. 

Our approximate A iteration procedure entails alternatively solving equation 26 for T"^^ 
and equation 21 for i^^^ . 

There are many possibilities for A*, and much effort has been invested in order to 
optimize its construction. Generally speaking, the choice of a particular A* is based on per- 
formance during numerical experimentation and varies from problem to problem. Following 
this heuristic approach we considered several different approximate operators and compared 
their convergence properties with the simple A iteration. Our test cases were the first 20 
time steps of two single grid collapse calculations with a) 13 x 13, b) 128 x 128 grid cells, 
and c) the results of time step 2900 of the 30 M0 case discussed in section 4. For these test 
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Table 1: TEMPERATURE ITERATIONS 



procedure 



a b c 

13x13 128x128 3*64x64 



A 



10.1 
10.9 
4.3 
4.35 



16.4 



57.8 




10.1 



> 100 



10.2 



Note. — Average over 20 test runs 

calculations we utilized the most efficient procedure for solving equation 21 for J"*'^ (dis- 
cussed below) . The average number of A iterations necessary for each approximate operator 
and for each test case is given in Table 1. 

Our approximate "core-wing" operator is constructed in analogy to the "core-wing" A 
operator of Scharmer (1981) which takes into account that lines are often optically thin in the 
wings and optically thick in the line's core. The analogy with a more or less monotonously 
decreasing (or increasing) continuum absorption coefficient is to assume a threshold value 
for which the dusty material becomes optically thin. We conducted several numerical ex- 
periments with different threshold values but found little improvement in comparison to the 
simple A iteration procedure. 

Defining A^'^^ = diag{A)~^ resulted in improved convergence behavior. Use of addi- 
tional accelerator terms as outlined by Auer (1987), however, did not always lead to further 
improvement. Our best convergence results were attained with the A^ operator constructed 
as follows. We first approximate Jij^k by using the diagonal operator A^'^^ in equation 24: 




(27) 



This approximation is used for the off-diagonal terms in equation 22 to obtain: 



Ji,j,k = (Oij^fc)"^ ji,j,k - Z)i ''^tij,kBi{Tij,k) 

~^'lj,k'^lj-l,ki'^) ~^ ^'lj,k'^l,j+l,ki'^) 
~^lj,k'^l,j,k-li'^) ~ ^ij,k'^l,j,k+li'^) 



(28) 
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Because the right hand side depends only on known quantities and the temperatures T this 
corresponds to the operator equation 

J, = A^.(T) . (29) 

Use of sometimes leads to oscillations during the iterations which we are able to damp 
by interjecting simple A iterations. 

Occasionally, the Newton- Raphson iteration procedure used for solving equation 26 did 

not converge. For these special cases we had to abandon our approximate A iterations in 
favor of the following more robust but much more CPU intensive procedure. If wc denote 
by Tj our (inaccurate) estimate of Tj then a temperature correction ATj can be determined 
from: 



where we have replaced the frequency integration by a weighted sum with the weights w^. 
By substituting By{Ti) + dB^/dT\rp^f,ATi for B^{Tj) on the right hand side of equation 20 
and replacing ATj by the expression given in equation 30 we find a modified FLD equation: 



V • (P.V J.) = J. + </ {-Ju - B,{f,) 



(31) 



which can be expressed in matrix form as 

A J. + 5^ Jm = (T) . (32) 

The iteration procedure now consists of alternatively solving equations 30 and 32 until con- 
vergence. Because this alternative procedure is more CPU-intensive than the approximate 
A iterations described above, it is only used when absolutely necessary. 



2.4-2. Solution of the partial differential equations 

The partial differential equations 20 or 31 have to be solved repeatedly for each hy- 
drodynamic time step and for each frequency. The coefficients V^,, and e^, depend on 
frequency and position and vary from time step to time step. The relative importance of 
individual terms can vary strongly from position to position and with time: i.e. both 

|V • V^VJ^\ > 



3u + Y.<:^\J'^-B.iTi)\ 
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and 

i 

occur. 

In addition to the "alternating direction implicit" procedure (ADI) we considered the 
multi-grid scheme "full approximation storage" (FAS), "successive over-relaxation" (SOR) 
and its special case by Gauss-Seidel (GS), the method of "quasi-minimal residues" (QMR), 
and bicgstab(2), an improvement on the "Bi-CGSTAB" algorithm ("Bi-CGSTAB" combines 
the advantages of GMRES and Bi-CG). More information on the ADI, FAS, SOR, and GS 
procedures can be found in Press et al. (1992). Our variant of the QMR procedure is 
described by Biicker & Sauren (1996) and the bicgstab(2) method is described in full detail 
by Sleijpen & Fokkema (1993). 

Our standard test problem was constructed from the density distribution calculated 
by Yorke, Bodenheimer, Laughlin (1995) for the case of a slowly rotating 10 Mq molecu- 
lar clump, 7074 years after the beginning of collapse (see their Fig. 4). The density was 

remapped onto a single 66 x 66 grid. Assuming a central luminosity of 100 and constant 
dust temperature (Tj = 20 K) as starting values, we constructed an initial model by iterating 
equation 20 until convergence using SOR, subsequently solving for the corrected tempera- 
tures using equation 18. The time required for this first iterative step was not included in 
our comparison. Each procedure was coded in fortranQO and optimized for vector processing; 
a single processor of a CRAY T90 was used for the comparison calculations, the results of 
which are shown in Figure 3 and Table 2 for ADI, QMR, bicgstab(2), and GS. 

Table 2: GPU TIME NEGESSARY FOR GONVERGENGE 



procedure 


1100 ^m 


435 nm 


ADI 


2.12 


> 600 


GS 


3.98 


14.30 


QMR 


2.17 


2.03 


bicgstab(2) 


1.15 


1.38 



Note. — Average over 10 tests. Calculations were run until the residual fell below 10 "^"ergs ^ cm Hz ^. 

We first note that the ADI procedure used by Sonnhalter et al. (1995) did not converge 
for several frequencies for our test problem; the residual in the central cell remained constant 
after a few iterations. This problem was alleviated by resorting to SOR; the most robust 



V • P^VJ^I < 
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bicgstab(2) 
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Q 
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Fig. 3. — Convergence properties of four versions of the partial differential equation solver 

for equation 20 arc shown at three wavelengths: 435 nm (high optical depth), 21 yum (optical 
depth of order unity), and 1100 /ini (small optical depth). The "residual" is the Euclidean 
norm of \AuJu — £i/(T)| over the entire grid. 



-20- 



variant used an over-relaxation parameter u — 1, reducing SOR to the GS procedure. GS 
iterations resulted in steadily (but slowly) decreasing residuals. QMR and bicgstab(2) re- 
quired the fewest number of iterations to reach convergence, but because of their complexity, 
it is not immediately apparent that these procedures would also require the least amount of 
CPU time. As shown in Table 2 bicgstab(2) did indeed prove to be the most efficient partial 
differential equation solver. It was used for our numerical simulations discussed in section 
4. Every variant of FAS we attempted diverged for our test problem and the method was 
quickly abandoned. 

3. Initial conditions 

In spite of recent observational and theoretical progress, the initial conditions for proto- 
stellar collapse are still poorly known. Whereas it is clear that the formation of massive stars 
require high masses, neither the average temperatures nor the sizes and density distributions 
of those molecular clumps on the brink of forming massive stars are especially well known 
(Stabler, Palla, & Ho 2000). 

Because we are investigating whether massive stars can form by accretion, comparable 
to our current understanding of low mass star formation rather than by coalescence of lower 
mass hydrostatic components or stellar coalescence (see e.g. Bonnell, Bate, & Zinnecker 
1998), we will adopt for our initial conditions a scaled- up version of the initial configuration 
expected for the formation of low mass stars (see e.g. Williams, Blitz, & McKee 2000 or 
Andre, Ward- Thompson, & Barsony 2000 for extensive reviews): clump sizes are a fraction 
of a pc, temperatures lie in the range 10-30 K (wc adopt Tq = 20 K), and the clumps are 
density-peaked towards the center (typically, p oc r^''\ where p ^ 1 —2). Because high mass 
star formation is an extremely rare event, we do not restrict ourselves to clump masses which 
are one or two Jeans masses only. We adopt a thermal to gravitational binding energy ratio of 
Et/\Eg\ — 0.05, corresponding to about 10 Jeans masses. However, because gravity always 
dominated thermal pressure forces in our domain of integration, we expect no significant 
evolutionary differences as long as Tq < 100 K (M > 2 Jeans masses). 

The initial configuration is summarized in Table 3. We begin with a rotating (Q = 
5 X 10~^^s~^) density configuration p oc r~^. This power law dependence corresponds to 
that expected in the central regions of magnetically supported, slowly contracting, slowly- 
rotating clouds (Mouschovias 1990; Tomisaka, Ikeuchi, & Nakamura 1990; Lizano & Shu 
1989; Crutcher et al. 1994). Probably the best measured pre-coUapse core, L1689B, has a 
somewhat flatter central region {R < 4000 AU) with p oc r""'^ or p oc r^^'^, depending on 
geometric and temperature assumptions (Andre, Ward-Thompson, & Motte 1996). Outside 
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Table 3: INITIAL CONDITIONS FOR COLLAPSE 



frequency-dependent 


F30 


F60 


F120 


grey cases 


G30 


G60 


G120 


Mass [Mq] 


30 


60 


120 


Radius [pc] 


0.05 


0.1 


0.2 


To [K] 


20 


20 


20 


Qo [10-1=^ s-i] 


5 


5 


5 


po [10-2°gcm-=^] 


1 


1 


1 


{uh) [lO^cm-^^j 


23 


5.8 


1.5 


Et/\Eg\ 


0.05 


0.05 


0.05 


Emi/\EG\ 


0.023 


0.094 


0.374 


ts [10^ yr] 


32 


65 


129 



Note. — The initial clump temperature Tq was also the outer temperature boundary condition Tout- 

4000 AU the density decreases as p oc r~^. By contrast, Motte, Andre, & Neri (1998) find 
p oc r^^ density profiles in 10 compact pre-coUapse cores in the p Oph region, which they 
were able to resolve down to ~140 AU. 

In spite of the fact that the radial dependence of density in pre-collapse cores is ob- 
servationally ill-constrained, especially for the high mass case, we have restricted ourselves 
to an initial p oc r^^ density configuration, corresponding to that of a "singular isothermal 
sphere" , which has been used for many years to model the formation of low mass stars. This 
allows direct comparision with the results of analogous collapse calculations. During the 
course of evolution, however, rotation, thermal forces at the outermost radius, and radiation 
forces quickly modify this initial density distribution. Moreover, in contrast to similarity 
solutions that rely on the assumption of a singular isothermal sphere our clumps are initially 
at rest. The mass accretion rates are thus not fixed, but rather a result of the radiation 
hydrodyanmic calculations. 

For the cases F120 and G120 the rotational velocity at the equator and Vmax exceeded 
the escape velocity. Thus, about 7.3 Mq of the original 120 Mq in the cloud was not bound 
initially. 
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4. Results 

4.1. Evolution of the Central (Proto-)Star 

We have implicitly assumed that all material flowing into the central zone is accreted 
onto a single object. Its mass is known from integrating the mass accretion rate over time; 
its luminosity and effective temperature are modeled according to the procedure described 
in section 2.1. In Fig. 4 we display the evolution of the total and accretion luminosities for 
the frequency-dependent "F" sequences. 

Considering the fact that the initial free-fall times double from case F30 to F60 and from 
F60 to F 120, the initial rapid rise of luminosity — due to the emission from an accretion shock 
surrounding the central star — appears similar for all three cases. This is not too surprising, 
because neither rotation nor radiation has an important influence during these early phases; 
the flow is dominated by gravity. After several thousand years, the energy released within 
the accretion shock front is not the dominant source of luminosity. Adding material onto the 
central star at the very high rates considered here causes it to bloat up to radii exceeding 
the deuterium burning hmit (see Fig. 5). Note that there is very little difference of the 
evolution within the HR diagram of these three cases. Indeed, except for a "shift" in time 
and the maximum mass attained in each simulation, the growth of mass within the central 
computational zone was remarkably similar (c.f. Fig. 6). The rate of mass accretion rises 
sharply within a few thousand years, reaches a maximum value ~ 2 x 10~^ yr~^. 

Comparing these tracks with those pubhshed by Behrend & Maeder (2001) by Meynet 
& Maeder (2000), we note that our tracks lie significantly higher in the HR diagram (larger 
radii, larger luminosity, but slightly lower Tcs for equivalent masses) and cross the main 
sequence at a somewhat higher luminosity by about 0.5 to 0.8 dex. This is not too surprising, 
because our luminosity includes 3/4 of the accretion luminosity (see equation 10) and our 
accretion rate varied in time according to the results of hydrodynamic calculations (see Fig. 
6 and Fig. 7). By contrast, Behrend & Maeder assumed a given accretion rate onto the star 
based on observations of outflows M* = max(10~^ M© yr~^, /Mout), where / was chosen to 
lie between 0.3 and 0.5. For outflow mass loss Mout the authors used the observed relation 
between the outflow mass rates and the stellar bolometric luminosities in ultracompact HII 
regions found by Churchwell (1998) and conflrmed by Henning et al. (2000). Meynet & 
Maeder used the formula = 10"^ M© yr'^ max(l, M^/M©)^ ^ 

The formulae assumed by Behrend & Maeder and by Meynet & Maeder result in higher 
accretion times (lower average accretion rates) and a shift of the phase of extremely high 
accretion rates to later evolutionary times, compared to what we flnd in our collapse calcu- 
lations. At the high mass end, however, the evolution of the star is similar; the tracks follow 
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Age [ 1000 yr ] 




Fig. 4. — Total luminosity {solid lines) and accretion luminosity {dashed lines) for the 
frequency-dependent "F" sequences 
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4.6 4.4 4.2 4.0 3.8 3.6 

log T^ff 

Fig. 5. — Evolution of the central (proto-)stars for the frequency-dependent "F" sequences 
in the HR diagram. The luminosity contribution Lace from the relaxation zone behind the 
accretion shock has been included. For comparison, the evolutionary tracks of non-accreting 
stars has been given (c.f. Fig. 1). 
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closely along the main sequence in spite of high mass accretion rates. 

I 1 1 1 1 1 1 r 




20 40 60 

Age [lOOOyr] 

Fig. 6. — Evolution of the central (proto-)stars' masses for the frequency-dependent "F" 
sequences. After an initial delay the mass accretion rates rise rapidly to comparable values 
(given by the slope of the curves) and fall off rapidly as radiative effects inhibit further 
accretion (see also Fig. 7). 

The growth of mass in the central computational zone was governed by the relative 
importance of centrifugal, radiative, and gravitative forces in the molecular clump. For 
the six cases calculated we found that the detailed treatment of radiation transfer strongly 
influenced the clump's evolution. In general, the accretion rate increased sharply after about 
one free-fall time and then decreases. To exemplify this we display the growth of central 
mass and the time dependence of the accretion rate for the frequency-dependent case F60 
and for the grey case G60 in Fig. 7. Assuming grey radiation transfer, in infall of material 
into the central regions of the molecular clump is strongly hampered after about 14 000 
years. The central star is a 20.7 M© main sequence star with a luminosity of 5.2 x 10^ Lq. 
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Assuming frequency-dependent radiation transfer, more than 60% more mass could accrete 
onto the central object. For case F60 the final mass was 33.6 Mq and the final main sequence 
luminosity was 2.2 x 10^ Lq. 



g20 - 
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Age [ 1000 yr ] 



1 r- 




\ 1 
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Fig. 7. — Evolution of central mass {left frame) and mass accretion rate {right frame) for the 
frequency-dependent case F60 {dashed line; shown previously in Fig. 6) and the grey case 
G60 {solid line). 



4.2. Evolution of Molecular Cloud Clumps 

In the following four figures the distributions of the gas density, of the temperatures of 
amorphous carbon and silicate grains and of the gas velocity are shown for selected cases 
(F60, G60, F30, and F120) at six selected times. The evolutionary age, total mass within 
the computational grid, total luminosity (including accretion luminosity), and mass of the 
central star which correspond to each frame are given in Table 4. 



4.2.1. CaseF60 

The evolution of the molecular clump considered in case F60 is shown in Fig. 8. Initially, 
the infall is nearly spherically symmetric. Even after 10 000 years (frame 8a) at a time when 
the core mass has grown to 13.4 Mq, there are only minor differences in the clump's structure 
along the equator compared to along the rotational axis, evidenced by slight flattening of 
density contours and elongation of temperature contours in the inner regions. After 20 000 
years (frame 8b) the non-sphericity has become more pronounced. In the polar regions 
directly above and below the central object (a 28.4 Mq star) the infall has been reversed by 
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Table 4: PARAMETERS OF EVOLVING CLUMPS 



case 


Age 


11 yf 


T 




Fig. 


[1000 yrj 


Fa /r 1 


[1000 LqJ 


[M©] 


rbU: oa 


iU 


DD.Z 


■1 1 


1 Q /I 

id. 4 


8b 


20 


73.2 


133 


28.4 


8c 


25 


72.9 


210 


33.0 


8d 


30 


62.2 


223 


33.6 


8e 


35 


47.0 


223 


33.6 


8f 


45 


39.9 


223 


33.6 


KjOU. ya 


iU 


00.4 


1 A 

14 


io. i 


9b 


25 


69.2 


52 


20.7 


9c 


35 


69.0 


52 


20.7 


9d 


90 


68.8 


52 


20.7 


9e 


100 


68.9 


52 


20.7 


9f 


110 


69.0 


52 


20.7 


roU: iUa 


1 n 
iO 


o/ .9 


z / 


ib. / 


10b 


15 


41 6 


85 


24.4 


10c 


19 


43.7 


143 


29.0 


lOd 


21 


44.6 


178 


31.2 


10c 


22 


44.5 


181 


31.4 


lOf 


24 


44.6 


184 


31.6 


F120:lla 


10 


125.3 


1.9 


8.1 


lib 


28 


134.0 


209 


32.9 


11c 


36 


134.6 


331 


38.3 


lid 


60 


121.3 


463 


42.9 


lie 


67 


109.1 


463 


42.9 


Uf 


76 


94.0 


463 


42.9 



Note. — 



Mgrid is the total mass within the computational grid, including the stellar mass M*. 
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radiative forces on the dust, The low density outflowing region is encased by shock fronts. 
The star continues to accrete material through the equator. After 25 000 years (frame 8c) 
the stellar mass has grown to 33.0 Mq and the cavity emptied by radiation has grown to 
±10^'' cm in the polar direction. Its expansion velocity has also grown (> 10 km s~^). 
Outside the cavity the infall of clump material has also been reversed; only in the equatorial 
plane is some infall possible. 

For the last three frames of Fig. 8 the central star has reached the main sequence at 
33.6 Mq and no longer accretes. Note that there is some indication of a flattened density 
structure (sec white contours in frame 8c), which persists for the remainder of the evolution. 
The highest expansion velocities (and lowest densities) occur towards the poles. Although ~ 
40 Mq of material was available within the computational grid at 25 000 years, radiative forces 
have effectively prevented significant accretion after this time. We stop the computations at 
an evolutionary age of 45 000 years (frame 8f). 

4.2.2. Case G60 

Although the evolution of case G60 initially proceeds in a similar fashion as case F60, 
notable differences are apparent after 25 000 years (compare frame 9b to frame 8c). The flow 
of material into the inner zone has effectively been stopped at 15 000 years, a time at which 
only 20.7 M© had been accreted (see Fig. 7). Contrary to case F60 there is no indication 
of the formation of a polar cavity, evacuated by radiative forces, even at the most advanced 
evolutionary time considered, 110000 years (frame 9f). Instead, the material flows onto a 
thin, disk-like structure, supported in the radial direction by both centrifugal and radiative 
forces. 



4.2.3. Case F30 

As for case F60 the initial collapse of case F30 is nearly spherically symmetric until an 
evacuated polar cavity is formed, encased in a system of expanding shock fronts. Here, the 
infalling material collides with the radiatively accelerated outflow. In the condensed cylin- 
drical shell bounded by the shock fronts (see frame lOe) the hard radiation from the central 
source is absorbed and reemitted at longer wavelengths. Because of this degrading of the 
stellar radiation "hardness" , the material outside the cylindrical shell is able to flow radially 
inwards more or less parallel to the equator for \z\ < 10^^ cm. Material continues to flow 
into the central zone via the equatorial plane, and the central stellar mass ultimately grows 
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Fig. 8. — Distribution of density (grey-scale and white contour lines)^ velocity (arrows)^ 
temperature of amorphous carbon grains {solid black contour lines) ^ and temperature of 
silicate grains {dotted contour lines) for case F60 (see Table 3) at evolutionary times as 
indicated in Table 4. 



-30- 





Fig. 9. — Distribution of density, velocity, and grain temperature for case G60 at evolutionary 
times as indicated in Table 4. Symbols and lines are as in Fig. 8. 
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Fig. 10. — Distribution of density, velocity, and grain temperatures for case F30 at evolu- 
tionary times as indicated in Table 4. Symbols and lines are as in Fig. 8. 
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Fig. 11. — Distribution of density, velocity, and grain temperatures for case F120 at evolu- 
tionary times as indicated in Table 4. Symbols and lines are as in Fig. 8. 
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to 31.6 Mq, more than originally present within the computational grid. The cylindrical 
shell depicted in frames lOd and lOe is a short-lived phenomenon, however. Within a few 
thousand years it collapses into a long, narrow, filamentary structure (frame lOf) containing 
about 13 MfT^. 



4.2.4. Case F120 

Again, as in cases F60 and F30 discussed above, the initial collapse is spherically sym- 
metric (frame 11a), followed by the formation of a polar cavity evacuated by radiative forces 
(frame lib), after a significant amount of material has accumulated within the central zone 
(M^ = 32.9 M0 at t = 28 000 yr.). The central star continues to accrete an additional 
10 Mq via an equatorial flow through a disk-like structure over the next 30 000 years albeit 
at an ever decreasing rate. This "disk" is short-lived, however. At an evolutionary age of 
60 000 years (frame lid) the accretion process has stopped and a cyhndrical shell first forms, 
contracts to a smaller cylindrical radius with a more focused polar outflow (frame lie), and 
then reexpands (frame llf). Whereas for earlier evolutionary times the gas density was 
higher in the equatorial regions than in the polar outflow regions, the flnal frame strongly 
resembles an expanding "cocoon" shell, punctured and elongated by a polar outflow. No 
disk-like structure is visible. 



4.3. Appearance of Molecular Clumps 

With the known density and equilibrium temperature distributions of each dust compo- 
nent calculated for each time step, we can perform ray-tracing radiation transfer calculations 
analogous to those of YB, who — by contrast — used the single dust temperature obtained in 
a grey radiation transfer code. From these ray-tracing calculations we can extract both SEDs 
(see Fig. 12 and 13) and isophote maps at selected wavelengths for any given evolutionary 
age. 

Because, however, our spatial resolution in the innermost regions is much worse than 
that of YB and furthermore, we have made no attempt to model the emission from a hypo- 
thetical disk within the central cell (see YB's "central zone disk model"), we will not be able 
to accurately model the emission from dust warmer than about 400 K. Even in the rather 
late formation stages considered here — the central stars have evolved to main sequence 
0-stars, very little near infrared and no optical/ultraviolet radiation escapes the remnant 
molecular cloud. It is reasonable to assume that a non-homogeneous distribution of material 
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Fig. 12. — SEDs at selected evolutionary ages of sequence F60. The four curves marked 
"a" correspond to the time shown in frame 8a, "b" to frame 8b, and "f" to frame 8f. The 
viewing angle 0° is pole-on, 90° is edge-on. 




Fig. 13. — SEDs at selected evolutionary ages of sequence F120. The four curves marked 
"a" correspond to the time shown in frame 11a, "b" to frame lib, and "f" to frame llf. 
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within the central zone (i.e. dumpiness or a flattened disk) would have allowed at least some 
of the hard photons from the central source to escape. Thus, the examples given here can 
only be illustrative of the basic method rather than an accurate model of the expected near 
infrared to ultraviolet flux. Unfortunately, our ignorance of the goings-on within the central 
zone preclude a more deflnitive treatment. 

5. Discussion and Conclusions 

Our improved frequency dependent radiation hydrodynamics code is able to track the in- 
fall of material within a molecular clump against radiative forces. We flnd that the "flashlight 
effect" first discussed by Yorke & Bodenheimer (1999), i.e. the non-isotropic distribution of 
radiative fiux that occurs when a circumstellar disk forms, is strongly compounded by the 
frequency dependent radiation transfer. The shortest wavelength radiation (which is also 
the most effective for radiative acceleration) is most strongly concentrated towards the polar 
directions, whereas the longer wavelength radiation (less effective radiative acceleration) is 
more or less isotropic. We conclude that massive stars can in principle be formed via accre- 
tion through a disk, in a manner analogous to the formation of lower mass stars. A powerful 
radiation-driven outflow in the polar directions and a "puffed-up" (thick) disk result from 
the high luminosity of the central source. 

We have developed a simplified model for following the evolution of accreting (proto-) 
stars, using existing tracks for non-accreting stars. With this model we have shown that in 
the case of massive star formation the energy released within the accretion shock front, the 
"accretion luminosity" , is not the dominant source of luminosity after a few thousand years 
of evolution. 

The accretion rate onto the central source is time- dependent. It rises sharply after one 
free-fall time to a maximum value and falls off gradually (in the frequency-dependent cases). 
This is in contrast to the expectations of Mcynct & Macdcr (2000) and Bchrend & Macdcr 
(2001), who have assumed mass accretion rates that increase monotonically in time up 
to a maximum value. 

We have also shown that the concept of "birthline", the equilibrium position of fully 
convective, deuterium-burning stars in the HR diagram with cosmic deuterium abundance, 
is — strictly speaking — unattainable for stars more massive than 1 Mq. Beginning with 
a protostar of a fraction of a solar mass and building up via accretion to 1 Mq and higher 
masses, it either accretes too rapidly (shifting the HR position to smaller radii) or it accretes 
too slowly (significant amounts of previously accreted deuterium are consumed) . For masses 
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M < 10 Mq, however, the contribution of the accretion luminosity may make the star appear 
to he on or above the birthhne. 

In this investigation we have not addressed the issues of the longevity of the circumstellar 
disk or the possible formation of a dense stellar cluster within our central computational zone 
rather than a single star. However, even without the assumption of ionizing radiation, we 
find that these disks are not long-lived phenomena. In the most massive cases the effects 
of radiative acceleration eventually disperse the remnant disks. Future studies will have to 
address the issues of ionization and the interactions of the disk with powerful stellar winds. 
The effects of nearby companions in a dense stellar cluster will also have to be considered in 
future work. 
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